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REMARKS 



The claims as filed were misnumbered beginning with claim 29 and the 
misnumbering is regretted. This has now been corrected. 

Claims 1-3, 5, 7-10, 13-22, 25-30, 36, 37, 40, 42,43,46-53, 56-58, 62-64, 66, 68, 
69, 72-79, 82, 83, 87-89, 91-94, 97-104, 106-108 and 1 12-117 are rejected under 35 
U.S.C. 102(e) as being anticipated by U.S. Patent No. 5,963,329 to Conrad et al. 
("Conrad"). Some of the rejected claims have been amended and the rejection is 
traversed in so far as it is applied tot he claims as amended. It is believed to be well 
established that for a reference to anticipate a claim, there must be identity of elements 
between the reference and the claim. 

Conrad fails this test as applied to claim 1 as amended. Conrad fails to teach or 
suggest detecting changes in polarization state of the beam that is diffracted by the 
diffracting structure, providing a set of polarization state data and arriving at the values of 
one or more parameters of the diffracting structure from the measured changes in 
polarization state and from an optimized estimation within a neighborhood of the set of 
polarization state data that is provided. The Examiner points to column 4, lines 8-15, 28- 
56 and column 8 lines 1-15 of Conrad as showing the above features claimed in claim 1. 
None of the sections referred to by the Examiner disclose the features in claim 1 . 
Column 4, lines 8-15 and 28-56 fail to teach anything even remotely resembling the 
above feature. In column 8, lines 1-15 Conrad describes a system for measuring 
parameters of a grating, where a polarizer 20 is used for providing polarized light column 
4, lines 8-15, 28-56 directed to the grating for measurement. The use of polarized light 
for the measurement by Conrad permits the grating profile to be calculated independently 
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from two normal polarizations: transverse magnetic (TM) that transverse electric (TE). 
As noted by Conrad, the grating profile calculated from these two different polarizations 
give almost identical line profiles. Column 10, lines 35-55. 

Since Conrad does not detect the polarization of light that is diffracted or reflected 
by the grating, it would be impossible for Conrad to detect the changes in polarization 
state caused by the diffraction/reflection. In contrast, such change in polarization state is 
measured in the invention of claim 1 . As clearly described in the embodiment on page 9, 
lines 12-22 of the specification of the present application, an analyzer 32 is used to detect 
the polarization state of the radiation diffracted by the diffracting structure, and the 
change in polarization state of such radiation caused by the diffraction can be measured 
by rotating either the polarizer 28 or the analyzer 32 when spectrometer 34 is detecting 
the diffracted radiation passed by the analyzer and from the diffracting structure 60 at a 
plurality of wavelengths. The changes in polarization state data at different wavelengths 
can then be derived by computer 40 from the intensities detected by the spectrometer. 
While this embodiment is useful in illustrating the feature claimed in claim 1, the scope 
of claim 1 is not limited to such embodiment. The polarization state of radiation 
diffracted by the diffracting structure is also measured in this embodiment, which permits 
the detection of change in polarization state of radiation caused by the diffraction. This is 
in contrast to Conrad, which does not measure the polarization state of radiation 
diffracted by the diffracting structure. From the above, it is evident that Conrad does not 
measure change in polarization state of radiation caused by diffraction or reflection by 
the grating. 
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From the above, it is evident that there is no identity of elements between claim 1 
and Conrad so that claim 1 is not anticipated by Conrad. As noted above, polarizer 20 is 
used by Conrad to enable two different independent calculations to be made, one in the 
TM mode and the other in TE mode, and not for the purpose for measuring the change in 
polarization state. There is therefore no reason or motivation to modify Conrad in order 
to measure change in polarization state caused by the diffraction/reflection by the 
diffracting structure. Claim 1 is therefore believed to allowable. 

As for claim 21, the Examiner appears to have ignored the limitation that "the 
wavelengths of the intensity or change in polarization state data in the one or more sets 
are chosen as a function of the properties of the one or more layers," where the 
diffracting structure measured is on or under such one or more layers when the 
measurement is carried out. This is not taught or suggested by Conrad. As explained in 
the embodiment of the invention on page 16, lines 1 1-23, polysilicon in the film stacks 
16a and 16c over or below the diffracting structure is substantially opaque at wavelengths 
from the ultraviolet range to about 380 nm. Therefore, if the intensity or change in 
polarization state data that are analyzed are only those for this wavelength range, the 
effect of variations in optical characteristics of the top or bottom film stacks on the 
measurements can be essentially ignored. This will simplify the calculations. For 
example, if the profile of grating 16b is measured at wavelengths from the ultraviolet 
range to about 380 nm, the measurements will not change despite variations in the optical 
characteristics of the film stacks across the grating. While this one example illustrates 
the functionality of the feature claimed in claim 21, the scope of claim 21 is not limited to 
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such example. Conrad fails to teach or suggest such feature and claim 21 is therefore 
believed to be allowable. 

The feature of claim 27 is illustrated by the embodiment of the description on 
page 17, lines 9-16. As clearly explained in such section and in Fig. 8, the diffraction 
intensity can change rapidly with a small change in wavelength at some wavelength 
values but changes rather slowly in other wavelength regions. To more accurately 
represent the shape of the curve, it is desirable to increase the density of the data samples 
within regions where the diffraction density changes rapidly with wavelength such as 
within regions 184 and 186. While such example illustrates the functionality of claim 27, 
the scope of claim 27 is not limited to such example. Conrad fails to teach or suggest any 
feature even remotely resembling such limitation. Claim 27 is therefore believed to be 
allowable. 

Claims 2 and 3 are believed to be allowable since they depend from allowable 
claim 1, and further on the ground of the limitations in these claims. Claim 2 adds to 
claim 1 the generation of a library of sets of changes of polarizing state data and 
comparing measured changes in polarization state to the sets of data in the library, to 
define the set of change in polarization state data that corresponds to the first set of values 
of the one or more parameters. This is not taught or suggested by Conrad in the sections 
referred to by the Examiner or any other section. 

Claim 5 is believed to be allowable since it depends from allowable claim 1 . It is 
further believed to be allowable on the ground that it adds the limitation that the first set 
of values of the one or more parameters is chosen as a function of sensitivity of the 
change in polarization state data to changes in the one or more parameters. Column 5, 
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lines 28-47 referred to by the Examiner in the Office action does not teach or suggest 
such feature at all. Such section in Conrad merely describes a slab model for predicting a 
percentage of energy that is diffracted and has nothing to do with choosing the value of 
parameters as a function of the sensitivity of the change in polarization state data of the 
diffracted radiation. 

Claims 7, 8, 25, 26, 29 and 30 are believed to be allowable since they depend 
from allowable claims. 

The Examiner rejected claim 9 in view of column 9, lines 9-45 and column 4, 
lines 30-50 of Conrad. We respectfully disagree. While Conrad mentions the calculation 
of eigenvalues in column 9, line 37, Conrad fails to teach or suggest storing the 
eigenvalues and using the stored eigenvalues for obtaining the value of one or more 
parameters of the diffracting structure. 

Models of the interaction between the beam of electromagnetic radiation and the 
grating are used for solving Maxwell's equations of such interaction. In the solution of 
the Maxwell's equations, the process decomposes the solution into a number of steps. An 
intermediate step in solving the equations is the derivation of the eigenvalues. The 
invention of claim 9 is based on the recognition that, for many diffracting structures, the 
initial steps for solving the Maxwell's equations are the same so that these steps do not 
need to be repeated every time but rather can be saved for future reference and reused, 
thereby saving time and effort in calculation. See for example the specification of the 
present application, page 13, line 21 through page 14, line 15. For this reason, the 
eigenvalues may be stored for future reference and can be reused to save time and 
calculation. This is not taught or suggested at all by Conrad. 
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From the above, it is clear that there is no identity of elements between claim 9 
and Conrad. There is also no reason or motivation for modifying Conrad to arrive at the 
features of claim 9. Claim 9 is thus also nonobvious over Conrad and is therefore 
believed to be allowable. 

Claim 10 is believed to be allowable since it depends from allowable claim 9 and 
because of the features of the claim. Claim 10 adds the limitation that the providing 
provides the model that accounts for the material(s) of the structure. This is not taught or 
suggested by Conrad. Column 5, lines 30-35 of Conrad referred to by the Examiner 
describes a slab model which does not take into account the type of material in the 
diffracting structure, unlike claim 10. 

The Examiner rejected claims 13-18 in view of the disclosure in column 6, lines 
6-30 of Conrad. We respectfully disagree. Column 6, lines 6-30 of Conrad describes the 
S profile used in models for predicting diffraction from the diffracting structure. The S 
profile is used as an approximation to the shape of the structure, and is thus radically 
different from and has nothing to do with the S-matrix in claims 13-18. The S-matrix is a 
mathematical representation that characterizes the reflection and transmission of radiation 
at an optical interface. For a more detailed explanation of the S-matrix please see the 
article enclosed by Lifeng Li, entitled "Formulation and comparison of two recursive 
matrix algorithms for modeling layered diffraction gratings," J. Opt. Soc. Am. A/Vol. 13, 
No. 5/May 1996. 

Claims 19 and 20 are believed to be allowable since they depend from allowable 
claim 9. 
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Claim 22 is believed to be allowable since it depends from allowable claim 21 . It 
is further believed to be allowable since it contains the limitation that the wavelengths of 
the intensity or change in polarization state data are chosen to reduce the influence of the 
properties of the one or more layers (above and/or below the diffracting structure) on 
deriving the value of the one or more parameters. When electromagnetic radiation is 
directed at the structure to perform measurements, the measurements will be affected by 
the effects of such layers on the measurements. Column 8, lines 15-35 referred to by the 
Examiner pertains to a process for normalizing the illumination intensity of the 
illumination source despite change in intensity of the source and has nothing to do with 
changing the wavelengths to reduce the influence of the properties of the film stacks on 
deriving the parameters. Claim 28 is believed to be allowable for substantially the same 
reasons as those above for claim 22. 

For reasons similar to those described above for claims 1, 9, 21, and 27, claims 
36, 42, 52 and 56 as amended are likewise believed to be allowable. These limitations 
are not taught in the sections of Conrad referred to by the Examiner (column 7, lines 60- 
column 8, line 14 and column 5, lines 38-50) or teachings in any other section of Conrad. 

Claim 58 pertains to an apparatus for finding a value related to one more 
parameters of a 3-dimentional diffracting grating. Conrad has failed to teach or suggest 
such feature. It is noted, for example, that all of the models illustrated in the figures or 
described by Conrad are two-dimensional models which cannot be used by themselves 
for deriving parameters of a three-dimensional diffracting structure. Furthermore, two- 
dimensional models cannot be extended automatically in deriving parameters of three- 
dimensional diffracting structures. Since three-dimensional diffracting structures include 
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more variables and parameters that would need to be taken into account, and the number 
of variables that need to be taken into account in two-dimensional is already quite high, it 
is not obvious how the larger number of three-dimensional diffracting structures can also 
be taken into account to extend the two-dimensional models to three dimensions. 

Based on the above, there is no identity of elements between claim 58 and Conrad 
and further that claim 58 is nonobvious over Conrad. 

As for claims 37, 53 and 57, it is noted above that column 8, lines 15-35 teaches 
nothing more than normalizing the intensity of the illumination beam despite variations in 
the intensity of light from the illumination source, which has nothing to do with the 
features in these claims. These three claims are also believed to be allowable for such 
reason and since they depend from allowable claims. 

Claim 40 is believed to be allowable since it depends from allowable claims and 
also on the grounds of limitations therein, that is, for the same reasons as those explained 
above for claims 5, 53 and 57. Column 5, lines 28-47 referred to by the Examiner fails to 
disclose any feature even remotely resembling that in claim 40. 

For the same reasons as those explained above for claim 10, claim 43 is also 
believed to be allowable. Column 4, lines 35-40 of Conrad fails to disclose any feature 
resembling such feature in claim 43. 

Claims 46-5 1 are believed to allowable for the same reasons as those discussed 
above for claims 13-18. 

Claims 62-64, 66, 68, 69, 72-79, 81-83, 87-89, 91-94, 97-104, 106-108 and 112- 
1 17 are believed to be allowable for the same reasons as those discussed above for other 
rejected claims. As for claims 112-117, these are believed to be patentable over Conrad, 
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even in view of Conrad's disclosure in column 4, lines 20-27 ', since these two claims 
require that a processing instrument alters its one or more processing parameters in 

response to the set of value(s). Conrad's disclosure in column 4, lines 20-27 is nothing 

v 

more than a very general statement that the values of the parameters such as line width 
obtained are "useful for monitoring or controlling a process." Nothing in such statement 
teaches or suggests that the values can be used for actually altering one or more 
processing parameters of the processing instrument such as a stepper or etcher. While 
not limited to the embodiment described in the specification of the present application, 
page 20, line 16 to page 21, line7, this feature of claims 1 12 and 116 does permit an 
integrated single tool to be constructed, which measures a diffracting structure created 
using the processing instrument, and uses feedback control to correct any errors in the 
processing instrument. 

Claims 4, 39, 65, and 90 are rejected under 35 U.S.C. 103(a) as being 
unpatentable over Conrad. The rejection is respectfully traversed. 

These four claims are believed to be allowable since they depend from allowable 
claims. They are further believed to be allowable since they add the limitation of 
performing non-linear regression which is not taught or suggested by Conrad. The 
examiner admits that Conrad fails to disclose such feature, but has failed to cite any 
factual support for concluding that this limitation would have been obvious in view of 
Conrad. The reason used in the rejection, namely that "this amounts to the use of a 
different or preferred method of data analysis," assuming arguendo to be true, still does 
not lead to the conclusion that this limitation, admitted by the examiner to be absent from 
Conrad, would have been obvious in view of Conrad. 
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We appreciate the Examiner's indication that claims 6, 1 1, 12, 23, 24, 38, 41, 44, 
45, 54, 55, 67, 70, 71, 80, 81, 95, 96, 98 and 105 would be allowable if rewritten in 



also believed to be allowable. 

Claims 1 18-139 have been added to more completely claim the invention; these 
claims are likewise believed to be patentable over Conrad and all other art of record. 

Claims 1-30, 36-58, 62-83, 87-108 and 1 12-139 are presently pending in the 
application. Reconsideration of the rejections is respectfully requested and an early 
indication of the allowability of all of the claims is earnestly solicited. 



independent form. This has not been done since the claims upon which they depend are 
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Formulation and comparison of two recursive 
matrix algorithms for 
modeling layered diffraction gratings 



Lifeng L| 



Lifeng Li 

Optical Sciences Center, University of Arizona, Tvcson, Arizona 85721 

Received July 20, 1995; revised manuscript received November 13, 1995; accepted December 4, 1995 
Two recursive and numerically stable matrix algorithms for modeling layered diffraction gratings, the S-matrix 
algorithm and the R-matrix algorithm, ara systematically presented in a form that is independent of the 
underlying grating models, geometries, and mountings. Many implementation variants of the algorithms 
are also presented. Their physical interpretations are given, and their numerical atabihties and effiaenciea 
are discussed in detail. The single most important criterion for achieving unconditional numerical atamlity 
•with both algorithms is to avoid the exponentially growing functions in every step of the mata* recursion. 
From the viewpoint of numerical efficiency, the S-matrix algorithm is generally preferred to the ii-matnx 
algorithm, but exceptional cases are noted. © 1996 Optical Society of America 



1. INTRODUCTION 

As research in the field of diffraction gratings advances 
and the range of grating applications widens, the struc- 
tures of gratings become more complicated than before. 
One of many new types of gratings that are finding more 
applications is layered gratings. For example, multi- 
layer thin films were deposited onto photoresist surface- 
relief gratings to make high-efficiency, all-dielectric 
reflection gratings, 1 and coating polycarbonate lamellar 
gratings with a layer of MgF* was proposed as a means 
of making broadband antirefiection structures. 2 Per- 
haps the most extreme cases of layered gratings are the 
Bragg-Fresnel gratings for use in x-ray spectroscopy 3 and 
the photonic band-gap materials* 4 On the other hand, 
in some grating models even a grating that consists of a 
single periodically corrugated surface is treated numeri- 
cally as a layered structure. In this paper the term 
layered gratings will be used broadly to refer to both 
physically and numerically layered periodic structures. 

All numerical methods for analyzing layered gratings 
face a common difficulty associated with the exponential 
functions of the spatial variable in the direction perpen- 
dicular to the grating plane- This difficulty is indicative 
of many problems of wave propagation and scattering in 
layered systems, and it is exacerbated by the fact that 
accurate numerical analysis of gratings usually requires 
a large number of eigenmodes. Recently this numerical 
difficulty has been overcome by many authors. 4 " 18 First, 
Pai and Awada 5 presented a Bremmer series method, 
based on the modal analysis with Fourier expansions, for 
dielectric gratings of arbitrary jprofile and groove depth 
in TB polarization. At the same time, DeSandre and 
Elson 6 presented an extinction-theorem analysis of 
diffraction anomalies in multilayer-coated shallow grat- 
ings by using the R-matrix propagation algorithm. 
Later, Unapplied the R-matrix algorithm to the classical 
modal method and enabled the latter to treat gratings 



of arbitrary profile, depth, and permittivity; Chateau 
and Hugonin 8 proposed an algorithm, with the coiipled- 
vrave method, to model surface relief and volume gratings 
made of lossless -and lossy dielectric materials. Montiel 
and Nevjere 9 applied the R-matrix algorithm to the dif- 
ferential method and thereby eliminated "the numerical 
instabilities that have plagued the differential theory in 
TM polarization during the past 20 years (Ref. 9, p, 3241). 
Recently Li 10 applied the R-matrix algorithm to the dif- 
ferential formalism of Chandezon et al (the C method) 
and thus removed a formerly existing limitation of the 
C method. The same goal was later achieved by Cotter 
e t aL lx using a scattering-matrix approach (S-matrix 
algorithm). The S-matrix algorithm was also used by 
Maystre 4 in an electromagnetic study of photonic band 
gaps by the integral method. Additionally, Ii w showed 
that under certain conditions the S-matrix algorithm 
(which, unfortunately, was referred to there as the 
R-matrix algorithm) and the Bremmer series algorithm 
are equivalent. Very recently . Moharam etal™ pre- 
sented another stable algorithm, which they call the 
enhanced transmission matrix approach. For references 
concerning the applications of the S-matrix end R-matrix 
algorithms to problems of wave propagation and scatter- 
ing outside the field of difi&action gratings, the reader 
may consult the reference list in Re£ 12. 

Now there exist many stable numerical algorithms 
and several variants of implementation, expressed with 
different terminologies and applied to different grating 
models. There are obvious similarities and subtle 
ences among these algorithms and their variants. Their 
advantages and disadvantages, as well as interrelation- 
ships, have not been addressed in the literature. The 
purpose of this paper is to provide a systematic and 
unified presentation of the S-matrix and R-matrix al- 
gorithms, independent of the underlying grating models 
(integral, cUfTerential, modal, etc) being used and the in- 
cidence conditions (TE, TM, or conical mount), and to 
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compare these two algorithms in terms of their physi- 
cal interpretations, numerical stabilities, and numerical 
efficiencies- Some results presented here have already 
appeared in the literature, but many intricate details 
are new. 

The algorithmic structure of the S-matrix and -K-matrix 
algorithms is recursive, and the matrix dimension in the 
recursion is independent of the number of layers. Mean- 
while, there exist stable and nonrecursive algorithms, for 
example, those in Kefs. 14 and 15. In these algorithms 
the field amplitudes in all layers are solved together from 
a large linear system of equations whose matrix dimen- 
sion is proportional to the number of the layers. The 
nonrecursive algorithms and the recursive algorithm of 
Mohararn ei a/., 13 which has a structure different from 
that of the S-matrix and E-rnatrix algorithms, are not 
considered in this paper. 

In what follows, first the framework is laid down, in 
Section 2, for the development in the subsequent sections 
by defining the notation and the basis functions. The 
S-matrix algorithm and the fl-matrix algorithm are. pre- 
sented in Sections 3 and 4, respectively. The presenta- 
tions axe arranged as parallel as possible for the two 
algorithms to bring out their similarities- Several vari- 
ants of the two algorithms are then given in Section 5, 
In Section 6 the two algorithms are compared in terms 
of their numerical stabilities and efficiencies. Finally, in 
Section 7 some remarks are made that are specific to the 
applications of the algorithms to several grating models. 

2. BACKGROUND FRAMEWORK 
A< Layer Abstraction 

Figure 1 depicts a multilayer surface-relief grating. We 
assume that the profiles of all medium interfaces have 
the same periodicity in the x direction and that they are 
invariant in the z direction. We say that two adjacent 
interfaces are separable if a line y = constant can be 
drawn between them without crossing «^r interface; 
otherwise, we say that they are nonseparable. Thus the 
bottom three interfaces in Fig. 1 are separable, and the 
top three are not 

The S-matrix and fl-matrix algorithms are applicable 
to all grating models, but here the discussion will be re- 
stricted to the classical modal, method, the C method, 
the coupled-wave method, and the differential method 
When it is not necessary to make the distinction, the first 
three methods will be referred to collectively as the modal 
methods, because they all rely on finding eigenmodes of 
Maxwell's equations. The classical modal method and 
the coupled-wave method approximate a continuous pro- 
file by a stack of lamellar gratings, as illustrated in the 
triangular grating in Pig. 1. This numerical approxima- 
tion effectively introduces a number of numerical layer 
interfaces. The differential method does not use the 
multilayer lamellar grating approximation, but for nu- 
merical purposes it decomposes a grating profile into thin 
horizontal slices, thus also creating numerical layer inter- 
faces. If two adjacent medium interfaces nave identical 
functional form and amplitude, the C method does not 
require any numerical layer interface; otherwise, one nu- 
merical interface may be needed between the two medium 
interfaces?^ 



Vol. 13, No. 5/May 1996/J. Opt. Soc. Am, A 1026 

We abstract a layered grating structure by a series of 
parallel straight lines, each representing a real or numeri- 
cal, straight or curved interface, depending on the profile 
of the medium interface and tne grating model being used 
[see Pig, 2(a)]. For example, suppose that for the layered 
grating shown in Fig. 1 we use the classical modal method 
to treat the rectangular profile, the same method with a 
three-layer approximation to treat the triangular profile, 
the differential method with a three-slice decomposition to 
treat the asymmetrical smooth profile, and the C method 
to treat the top three profiles. Then, in Fig. 2(a), n = 15, 
2 for the rectangular profile, 4 for the triangular profile, 
4 for the asymmetrical smooth profile, and 6 for the top 
three profiles. The permittivities in Fig. 2(a) may be ei- 
ther constants or periodic functions of x, depending on the 
spatial region and the grating model. Media 0 and n + 1 
are two semi-infinite homogeneous media. The dashed 
line in medium 0 is a numerical interface. It can be ar- 
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Fig, X. General layered grating. All periodic medium inter- 
faces share a common period, but otherwise they are arbitrary. 
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Fie. 2. Abstract layered grating structure, where the horizontal 
£L ac&l^aterlal interfaeea or nujner^ 
interfaces. The fields in each layer can be represented wther 
Ufa* a superposition ef upward, and da^ward-prop^ti^ 
and decaying waves or <W as a superposition of two sets of 
orthogonally polarized eigenmodes. 
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.trarily close to interface 0- The thickness of layer p 
ill be denoted by h p - 

I Basis Functions , , 

When apodal method is used to analyze layer p in 
Fie 2(a). the fields there are represented by superposi- 
tions of the eigenmodes. We assume that the eigenmodes 
in all layers share a common Floquet exponent, which is 
determined by the incident plane wave. The eigenvalue 
spectrum a™, having elements jSP, can be partitioned 
such that <r<» > = U o- 1 "" , where 

(1) 

In general, for any numerical truncation, a^" and «r lp) " 
have the Bame number of elements. The dependence of 
the eigenmodes ony. i.e., the y-dependent basis functions, 
is given by expiry], where J&f* e <r"»=. (Here y 
should be replaced by u if the C method is used; how- 
ever, for simplicity we will ignore this minor difference.) 
Thus we call an eigenmode corresponding to Am .an up 
wave, and that corresponding to A*™" 5 " a down wave. In 
particular, in the two semi-infinite regions and the ho- 
mogeneous regions between the separable medium inter- 
faces, the eigenmodes are simply the Rayleigh modes. In 
Fig. 2(a) the upward and downward arrows schematically 
represent the up wave and down waves, and the boldface 
letters u and d denote the column vectors whose elements 
are the wave amplitudes. Once the eigenmodes are de- 
termined everywhere, the grating problem reduces to a 
problem of determining the mode amplitudes. 

Alternatively, we can choose cos(Am y) and sin(Am y) 
as the basis functions, which is always possible because 
A 1 /'" = -A ( m p,+ with the classical modal method and the 
coupled-wave method and with the C method when the 
grating profile is symmetrical. We use U and V to denote 
the amplitudes of the modes that use this basis function 
set The physical meaning of these amplitudes is clear: 
for example, if U is proportional to the z component of the 
electric field, then V is proportional to the x components 
of the magnetic field, as«hematically shown in Fig. 2(b). 
From a mathematical point of view, the derivative of a 
U mode is a Vmode, and vice versa. For the.modal meth- 
ods both the exponential (or u-d)"basis functions and the 
trigonometrical (or U-V) basis functions can be used. 

In the differential method, one does not seek the eigen- 
solutions of Maxwell's equations. Instead, one numeri- 
cally integrates V from one interface to another, where U 
is a column vector whose elements are the Fourier expan- 
sion coefficients of the z component of the electromagnetic 
fields The numerical integration procedure gives the 
values of V and V = dI7/dy as functions of y. Clearly, 
the U-V basis functions described here correspond to 
the U-V basis functions described in the preceding para- 
graph. Thus the schematic diagram in Fig. 2(b) applies 
to the differential method as welL It is possible to have 
a set of u-d basis functions for the differential method 
if suitable linear combinations are made. 8 It is pnpor- 
tant to realize that, although the basis functions for the 
differential method do not have an explicit y dependence, 
their y dependence is asymptotically the same as that of 
the basis functions in the modal methods. 
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C. Boundary Conditions 

In this subsection we affix the equation numbers of all 
equations that apply only to the u-d basis functions 
with a letter a, and those that apply only to the U-V 
basis functions with a letter b. The Bame convention will 
be used for the formulas of the S-raatrix and R-matrbc 
algorithms in Sections 3, 4, and 5. 

In the modal methods, when the boundary conditions 
are matched along interface p, we generally get an equa- 
tion of form 

wlp ^^yp + « 1 « w4 U ! P ! ( * " (2a) 
W \d<>+Hy P + 0)\ W L d<P> U>- 0 >J 

where W lp) and W lp+V are square matrices. Further- 
more, by virtue of the modal fields, 



(3a) 



where 



rul'Hy, - 0)"| _ JV-to-i + 0)1 



.<,,., [exp^M 9, I- (4 

* 0 expO^-Uj 



a) 



and the exponential functions represent diagonal matrices 
(henceforth, all quantities with a subscript m represent 
diagonal matrices). Thus we obtain a recursive relation 
for the field amplitudes 

U^Hy P +o)l^ f J^Hy P -i^o)\ 

|_d<* + »(y p + 0jj Ld<«(y P -i + 0)J 



where 



with 



ftp) m Ty^+M-W' . 



(6) 



(7) 



The matrices and P' 3 can be fittingly called interface 
and layer t matrices, respectively. Note that |«» is of 
order 0(1) » iViU 

If the U-Vbasis functions are used to match the Douna- 
ary conditions, we have, correspondingly, 



W lP+U 



( ru<> + «(y p + 0)l w J Ui y>>- Q) l (2b) 

r cos^A,) 
L-^"- 1 sin(Aif 

[W+Hyp + 0)1 p p) rP , "U- 1 + 0)1 (5b) 

[Vl' +1 >(y, + 0)j 1 L V ° ,) ^ + 0) J 

In Eq- (4b), tj#' is a constant independent of h P , and for 
simplicity the plus sign has been dropped from tne su- 
perscript of eigenvalue a1? )+ . Also, the same notation 



,Jr»- 1 ■nCAlf 1 *,) cosfA^'M J 
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W 9 4> % and t is used with the two different basis function 
sets. This, however, is not a problem because the context 
will tell to which basis the matrices axe referring- From 
now on, the amplitude vectors without an explicit argu- 
ment stand for the vectors that are evaluated at the lower 
bound of the layer. For example, = W^jfc-i + 0). 

At this point, it is most natural and mathematically 
simplest to proceed with solving the grating problem by 
the so-called ^matrix algorithm, which is obtained by re- 
peated use of Eq. (5a) or Eq. (5b). However, it is well 
known that the r-matrix algorithm, with either of the ba- 
sis function sets, is numerically unstable when the total 
layer thickness of the grating structure and the matrix 
dimension are large. 7 This numerical instability is gen- 
erally attributed to the presence of the growing expo- 
nential functions in the algorithm. Fundamentally, the 
cause of instability is a classic one: loss of significant 
digits when one is computing a small number by subtract- 
ing two large numbers by a computer of finite precision. 
Symbolically, it is a case of » - x » 0(1). It should be 
emphasized that the numerical instability of the T-matrix 
algorithm cannot be eased or removed by simply reduc- 
ing the individual layer thicknesses without lowering the 
total thickness, because the T-matrix algorithm accumu- 
lates the magnitudes of the exponential functions as the 
layer r matrices are multiplied together. 

3. S-MATRK ALGORITHM 

The S-matrix algorithm uses the exponential basis func- 
tions. For any 0 £ p ^ n, it seeks a stack S matrix,. S ip \ 
that links the waves in layer p + 1 and medium 0 in this 
way: 

d^ +11 J' 

Before moving on, it is important to describe the physical 
meaning of the S matrix. For this purpose we rewrite 
Si p) in a two-by-two block form: 



[ d'°» J L d " 



(8a) 



d<°> J 



«,</>> 

t<Ld 



d^ +1) J' 



(9a) 



The significance of the subscripts u and d becomes evident 
once the reader mentally carries out the matrix-vector 
multiplication on the right-hand side of the equation. As 
a rule in this paper, the use of subscripts u and d is 
an automatic indication that the submatrix belongs to a 
matrix in the S-matrix algorithm. The choice of letters 
R and T, instead of S, makes the physical meanings of the 
four submatrices of S ip} self-explanatory. For example, 
Tiu and R { /d are the transmission matrix and reflection 
matrix that give the upward wave amplitudes in layer 
p + 1 resulting from the transmission of the upward 
incident waves in medium 0 and from the reflection of 
the incident downward waves in layer p + 1, respectively, 
by the whole stack below layer p + 1- Alternatively, the 
first p layers of the layered grating can be viewed as a 
linear four-terminal network. Matrix S M operates on 
the two sets of inputs to generate the two sets of outputs, 
as shown schematically in Fig. 3(a), 



\ \ 












o 




R (P) 





CP* f& 
(a) 



Fig. 3. Schematic diagrams for alternative interpretations of 
(a) the S matrix and (b) the R matrix. In Fig. 3(a), 1 and 0 
stand for inputs to and outputs from the system represented 
by the square box. In Fig. 3(b), i and v stand for currents 
and voltages at the terminals of the circuit represented by the 
square box. 

To link the waves in two adjacent layers, we can define 
an interface s matrix, s (p, p and a layer s matrix, as 
follows: 



ru<'-»(*+0n , p) r u^-o i 

Ld^(yp-O) J Ld^(y P + 0)J 



or 



U<" J 



S_]L*'' +1, J" 



(10a) 
(11a) 

(12a) 



The physical interpretations of s^' and a'' 1 are similar 
to that of S ( p\ Note that and because of the 
notation of their subscripts, cannot be confused with the 
t matrix defined in Section 2. The layer s matrix is re- 
lated to the interface s matrix by 



sM 



(13a) 



(14a) 



I" 1 ° lJe*9VM% 

[0 exp(-fAi?'"/ip)J 

and the interface s matrix ia in turn related to the inter' 
face t matrix by 

Note that all four entries in Eq. (Ha) contain the inverse 
of submatrix t<& ] . For this reason we call fj£* the pivotal 
submatrix. 

From Eos. (8a) and (11a), the set of recursion formulas 
for the stack S matrix are 

r^-^i-^r^T 1 ^ 
itf-itf+wifi-li-^airT 1 ^ ' 
& = *ir + lirw [i - Bir^'Vis-" 

Tiff - ^"[i - ijMrVtf; (15a) 
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The S-matrix recursion can be initialized by setting 

^"-[i a6a) 

or, equivalently, by setting S tOT « s i0 K 

The form of the factors enclosed by the square brack- 
ets in Eq. (15a) is such that the inverse matrices can be 
readily expanded* at least formally, into a geometrical se- 
ries in terms of the product of the two reflection matrices. 
This fact naturally gives rise to the multiple-reflection in- 
terpretation of the S-matrix recursion formulas. 1 * (It Is 
quite unfortunate that in Ref. 12 the S-matrix algorithm 
was incorrectly called the i?-matrix algorithm.) Because 
of the elegant form of the inverse matrices, Eqs. (16a) will 
be called the normalized S-matrix recursion formulas. 

Equations (8a)-(16a) constitute the basic ingredients 
of the S-matrix algorithm. In most grating problems the 
quantities of interest are the field amplitudes leaving the 
grating structure in the two outer media, i.e., u (n+1) and 
d <0j - They are simply given by 

r^iJis *3ir i (Ha) 

In particular, if there are no incident waves in medium 

ir = J$d< ft+l) , 

d (0) = T w d cnM, (i 8a) 

The numerical stability of the S-matrix algorithm 
is rooted in the construction of the layer. s matrix. 
The problem-causing, growing exponential function, 
exp(iA*f*~A,), that was originally in the t matrix, is 
now inverted in Eq. (13a). Since s (p * is of order 0(1), so 
is s**>K Furthermore, the submatrices of s ip) appear in 
the recursion formulas only as additive or multiplicative 
terms. Thus the S matrices remain of order 0(1), and 
the numerical stability of the algorithm is ensured. 

4. B-MATRIX ALGORITHM 
The E-matrix algorithm uses the trigonometrical basis 
functions. * For any 0 ^ p s n, it seeks a stack R matrix 
R { p\ that links the fields in layer p + 1 and medium 0 
in this way: 

The R matrix can be physically interpreted as field 
impedance or field admittance (the ratio of the tangential 
component of the E field to the tangential component of 
the H field or its inverse). For example, in TE polariza- 
don, because U and V correspond to the E and H fields, 
respectively, R xp) plays the role of field impedance. An 
alternative interpretation of Eq. (8b) can be made, with 
the aid of Pig. 3(b), in terms of currents and voltages In 
an electrical circuit. Here, if U is identified with the 
voltages and V with the currents, or vice versa, then R Lp) 
is the electrical impedance or admittance. The concept 
of impedance has been used previously in modeling grat- 
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ings that contain multiple planar interfaces but only one 
periodically modulated interface. 19 

To link the fields in two adjacent layers, we can define 
an interface r matrix, r tp K and a layer r matrix, r* p) , as 
follows: 

The layer r matrix is related to the interface r matrix by 
X{W v tfcs<{Wh p ], (13b) 

where 

For a proof of Eq, (13b), see Appendix A- Here we have 
excluded the possibility that accidentally )!mh p ». In, 
where Z is an integer.. The interface r matrix is in turn 
related to the interface t matrix by 

w j&4ir ^-W4 p) l. (14b) 

From Eqs. (13b) and (14b) it is clear that r ip) is of order 
0(1). ^Alternatively, we can use the layer f matrix to 
express the layer r matrix: 

(p) _ hi tzi tit -tii *2i in L (i4V) 

L <i? hl -ff^tf J 

In Eqs. (14b) and (14V) t[i and t£ are the oivotal sub- 
matrices. In some cases it is possible that t 2 i — 0, but 
then fii ) * 0. In fact, such a case can be utilized bene- 
ficially (see Appendix B). 

From Eqs. (8b) and (lib) the set of recursion formulas 
for the stack R matrix are 

Rl?~?\?-^Z<P>r[t 
R\? = r$zWR[r". 
R$ - -Btr"M**tit* 

r£ = asr 11 + Rir v z {p) R\£~", < l5b) 

where 

The J?-matrix recursion can be initialized by setting 

jjW^fW^ (16b) 

Unlike their counterparts in the S-matrix algorithm, 



9oota 



T9T9 $2,8 80* XVd CZ:ZT frOOZ/ZT/CO 



Ufeng Li 

Eqs. (15b) do not readily subject themselves to an intui- 
tive physical interpretation. Nonetheless, to preserve 
the formal symmetry between the two algorithms, we 
shall call Eqs. (15b) the normalized R -matrix recursion 
formulas. 

Starting with Eq, (16b), repeated use of Eqs. (15b) until 
p <= n leads to 



(17b) 



Suppose that = u«» + &*\ V« = u ,0) - d", t/ l * +1 ' = 
u 0it1i + d* ,l+1 > and V tn+1) = u <ll+11 - d< n+1 \ which is al- 
ways possible by definition. Then from Eq. (17b) we get 
the linear system that determines the out-going diffrac- 
tion amplitudes in the top and bottom media: 



L-j# i * j* 1 JL J 



■t 



■"12 



The jR-matrix algorithm is also immune to the numeri- 
cal difficulties associated with growing exponential func- 
tions- This is because the submatrices of r< pl are all of 
order 0(1), and they appear in the recursion formulas 
Eqs. (15b) and (15V) only as additive and multiplicative 
terms. The former fact is evident when Eq- (13b) is used 
to construct the layer r matrices. It is not so obvious 
if Eqs. (14b') are used, however; in fact, in this case the 
.R-matrix algorithm is only conditionally stable. 

If can be shown that Eq. (14b 7 ) and Eqs. (13b) ^are 
algebraically equivalent Therefore the submatrix * 
as given by Eq. (141/), should be mathematically pro- 
portional to csc(a£ 5 M, which tends to 0 as 771—°^ and 
v On the other hand, the first term of hi in 

]£. (14b') is - tW dn[AL%] + *!? cc*[^%], 
which tends to Thus the second term of ? t Z must also 
tend to * as m -» and h p -* »• Clearly this mathe- 
matical arrangement presents a serious numerical prob- 
lem. When the absolute values of the imaginary parts 
of I&f } h p are large, the numerically calculated matrix ele^ 
rnents of rtf by Eq. (14b 1 ) may not be small, as a result 
of round-off errors. Let A<* } be the maximum of the ab- 
solute values of the imaginary parts of all eigenvaluesfor 
a given matrix truncation. As a rule of thumb, when 
expLA^Ap] ~ 10 1 *, the numerical problem described 
above begins to arise (double precision is assumed here). 
To avoid the problem one has to choose the layer thick- 
ness so that exp[A"%] « 10 lS . Therefore thefi-matnx 
algorithm is conditionally stable when Eq. (14b') is used. 
Fortunately, unlike the T-matrix algorithm, here the mag- 
nitudes of the exponential functions do not- accumulate. 
Therefore lowering the individual layer thickness is an 
effective remedy for the numerical instability caused by 
the use of Eq. (14b')- 

5, VARIANTS OF IMPLEMENTATION 

A, Variation in Matrix Manipulation 

InSections 3and4theS-matrixandi?-inatrixalgoriOime 

were systematically presented. The presentation took 
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three steps: the definitions and derivations of the layer 
t matrices, the layer s (or r) matrices, and the stack S 
(or -R) matrices. Although from a theoretical point of 
view the introduction of the layer t matrices and the 
layer s (or r) matrices has made the presentation system- 
atic, from a practical point of view the use of one of the 
two kinds of layer matrices can be eliminated, as demon- 
strated below. 

The S-matrix recursion can be accomplished by use 
of the t matrices directly, without the layer 5 matrices. 
From Eqs. (6a), (6), and (9a) we can easily derive a set of 
nonnormalized S-matrix recursion formulas by using the 
interface t matrix: 

rtf - T^r V"- 1 ^* + n< " ) ]" 1 ' (19a) 



where 



(19a') 



and 4> { - fi) are the two diagonal submatrices in Eq. (4a). Of 
course, Eqs. (19a) and (16a) are algebraically equivalent. 
Note that the above equations have been written in terms 
of the interface t matrices, instead of the layer t matrices, 
and the appearance of 4> { I* has been arranged properly 
so that there are no exponentially growing functions in 
the formulas. This measure avoids possible numerical 
overflow and ensures numerical stability of the S-matrix 
algorithm. 

If we set 4>i p) = 1 in Eqs. (19a) and (19a') and replace 
all interface t submatrices by layer t submatrices, we ob- 
tain the nonnormalized S-matrix recursion formulas by 
using the layer t matrix: 

The use of this set of recursion formulas should be avoided 
whenever possible, because the matrix to be inverted is a 
sum of an exponentially growing matrix and an exponen- 
tially decaying matrix. . 

The ^-matrix recursion can also be accomplished with 
u6e of the t matrices directly, without the layer r matrices. 
Prom Eqs. (5b), (6), and (8b) we can easily derive a set of 
nonnormalized fl-matrix recursion formulas by using the 
layer t matrix: 

Since Eqs. (20b) are algebraically equivalent to Eqs. (15b), 
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*he R matrices obtained this way are mathematically of 
,rder 0(1). Numerically, however, devastating round-off 
errors could occur if the numerical layer thicknesses are 
set too high. The reason is the same as the one given at 
the end of Section 4 for the possible numerical instability 
resulting from the use of Bq. <l4b'). Specifically, the 
expression of*!? in Eqs. (20b) is of type « . - - - 0(1). 
Thus Eqs. (20b) also give a conditional stable unplemen- 
tation of the ^-matrix algorithm. The unconditional 
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B. Variation in Recursion Order 
The S-matrix and R -matrix recursions do not have to be 
performed in the order indicated in Sections 3 and 4, In 
other: words, one does not have to start the calculation 
Trom medium 0 and work step by step up to medium 
n + 1. For the normalized recursions this point can be 
best illustrated by the use of Redheffer's start product. 20 
Let a, 6, and c be 22V x 2N matriceE. Then the star 
product of a and b, in the S-matrix algorithm, is defined as 



*u\ la*+*«bi.it-a«b d J-*a Wi 0«(1 - b du a ltd ) b «* J 



stable, nonnonnalized R-matrix recursion formulas ob- 
tained by using the interface t matrix are given in 
Appendix C, , - , 

We recall that in Subsection 2.C we formally derived 
the layer t matrix from the boundary equation, Eq. (2a) or 
Eq, (2b). In fact, in at least two important cases Eg. (5a) 
or Eq- (5b) is obtained without the aid of Eq* (2a) or 
Eq. (2b). The firet case is the classical modal method in 
which matrix t ip) is obtained directly by projecting the 
functional boundary equations onto a natural basis func- 
tion set, 7 and the second case is the differential method 
in which matrix is obtained from a numerical inte- 
gration procedure. 3 In these cases the use of the non- 
normalized recursion formulas may be beneficial because 
they require feweT matrix operations. 

In the C method and the coupled-wave method the 
boundary equations (2a) or (2b) are an integral step of 
the numerical treatment In this case we can bypass the 
t matrix and derive the s (or r) matrix direcdy from 
the boundary equations. Writing the two W matrices m 
Eq. (2a) in a two-by-two form, and rearranging the terms ■ 
slightly, we have 



[war 1 -*sr 



J[dtP + «(y P +0)J 



(21a) 



Therefore from Eq. (10a), 



(22a) 



The layer s matrix, then follows immediately from 
Eq. (13a). Similarly the interface r matrix, rj* p can bo 
derived directly from boundary equation (2b); ie., 



(22b) 



The layer r matrix, f then follows immediately from 
Eq. (13b). 



where a^, b udy etc., are N X N submatrices. Similarly, 
we define the star product of a and b in the 22-matrix 
algorithm as 



fan 612 1 

ia Z i 022] \_bii b 22 J 

_ \b n - 612(^22 - a u )- l 6 2 i 012(022 " an)' x au 1 
L (23b) 

It can be shown that the star multiplications are associa- 
tive, i-e. t that 



a * (b * c) = (a * b) * c , 
a * (6 * c) = (a * b) * c , 



(24a) 
(24b) 



In the remainder of this section, for simplicity, I will 
mention only the 5-matrix recursion. The results for the 
H-matrix recursion can be obtained by obvious substitu- 



tions. 



In terms of star products, the S-matrix algorithm can 
be succinctly expressed as 



;(n) 



(25a) 



However, because of the associativity of the star multipli- 
cation, the product can be regrouped as follows: 



gin) =, 5 t<w * {..-. * [§to- 2 > * (£ tn *" 1) * * w )]- 



(26a) 



Equation (26a) describes the S-matrix recursion in the 
reverse order, starting from medium n + 1 and moving 
downward to medium 0. , , 

The associativity of the S-matrix recursion can be used, 
advantageously to save computation effort or to increase 
computation speed- Suppose that a large number of cal- 
culations are to be carried out for a grating with a vary- 
ing parameter that affects only thejth layer, for fixed j. 
Then the S-matrix recursion can be performed like this; 

S ("> [s» * ... • s"- u ] * * [S iJ ' l) * — • (27a) 

where the two recursions inside the square brackets are 
performed just once and then are used repeatedly to form 
a star product with the changing I ' * If the grabng cal- 
culation is done on a computer capable of parallel pro- 



900(g] 



H03M3X Ym 



T9T9 $2,9 $0t XVd f00Z/?T/C0 



Lifeng U 

ceasing, then the £ matrices can be grouped pairwise to 
increase the computation speed. 



6. COMPARISONS 

Having systematica^ addressed . the S-matrix and 
i?-matrix algorithms and their variants, we are now 
ready to make some comparisons- As mentioned in Sec- 
tions 3 and 4, the S matrices are related to the physical 
concept of reflection and transmission, and the R ma- 
trices are related to the physical concept of impedance 
and admittance. Furthermore, the normalized S-matrix 
recursion formulas can be readily interpreted in terms of 
multiple reflections, but the >R-matrix recursion formulas 
and the nonnormalized S-matrix recursion formulas are 
not easily interpreted in physical terms. In what follows 
we compare the numerical stabilities and efficiencies of 
the two algorithms. 

A- Numerical Stabilities 

Both the S-matrix and fl-matrix algorithms are inher- 
ently stable because the S matrices and the R matrices 
are both mathematically of order 0(1). However, there 
are subtle numerical differences between them. In gen- 
eral the S-matrix algorithm is much easier to work with 
than the R matrix algorithm. 

The implementation of the S-matrix algorithm is mostly 
worry free (but see Subsection 7.D), thanks to its use of 
the exponential basis functions. All s and S matrices are 
of order 0(1), not only mathematically but also numer- 
ically [we assume that Eqs, (20a) are not used). Thus 
there is no limit in layer thickness. The possibility of 
numerical overflow associated with the exponential func- 
tions is eliminated because only the decreasing exponen- 
tial functions are evaluated. Underflow can happen, but 
it is not a problem for most compilers. Additionally, the 
occurrence of Am *h p = Itr is not a problem at alL 

The implementation of thefl-matrix algorithm requires 
special treatment- Although all r and R matrices are of 
order 0(1) mathematically, they may not be so numeri- 
cally. When the factorization of the layer t matrix into 
the product of the interface t matrix and the diagonal ma- 
trix 4> is possible, one should use Eqs. (Cl) and (C2) below 
to perform the nonnormalized recursion or use Eqs. (13b) 
to compute if the normalized recursion is to be used. 
When the factorization is impossible, as is the case for 
the differential method, the numerical layer thicknesses 
should be kept sufficiently low that the computation of 
f\i by Eq. (14b') or of R[ p) by Eq, (20b) will not suffer 
significant loss of accuracy. There is also a minor tech- 
nical complication. Functions cot and esc that admit a 
complex argument are not intrinsic functions in most com- 
pilers. Thus the programmer has to write cot and esc as 
user-defined functions, using either the sine and cosine 
functions or the exponential functions. In doing so, care 
has to be taken to avoid overflow. In principle, the ac- 
cidental occurrence of Am'/i, = In is a problem for the 
^matrix algorithm. In practice, it is highly unlikely 
that the equality holds for / * 0 with high precision; 
therefore the singularity never poses serious numerical 
problems. Of course, one should judicially avoid hp.** 0, 
which is an uninteresting case, anyway. 
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B. Numerical Efficiencies 

In this subsection we consider the numerical efficiencies 
of different variants of the S-matrix and A-matrix algo- 
rithms. More specifically, we estimate the number of al- 
gebraic operations that each variant takes to compute the 
outgoing diffraction amplitudes u (n ' Pl) and d^, assuming 
that the W matrices in the boundary equations have been 
obtained. 

As is evident from the presentations in Sections 3 and 
4, after the S-matrix recursion is completed. u tn *« and 
d l0) are readily given in a solved form. However, with the 
fi-matrix algorithm the completion of thcR-matrix recur- 
sion only gives a system of linear equations that has yet 
to be solved to yield u tfl+l1 and d m . This initial compari- 
son is already in favor of the S-matrix algorithm 

We now take a closer look at the structure of the ma- 
trix recursion formulas. We say that a subset of the four 
submatrices of S ip) is a closed set with respect to the 
S-matrix recursion if every element of the set is deter- 
mined by the elements in the same set. Thus S ipi has 
four closed proper subsets: 



{Rud)t 



(28a) 



If there are incident plane waves in both media 0 and 
7i + 1, then from Eq. (17a) the S-matrix recursion of 
all four submatrices has to be performed. Suppose that 
u<°> = 0; then only the recursion of {R$, T$) is neces- 
sary if both u* n * 1> and d<* are needed, and only the recur- 
sion of R l Jd is necessary if only u« n+1) is needed. We call 
the recursions above the full, half, and quarter S-matrix 
recursions, respectively. 
Similarly, R ip) also has four closed proper subsets: 



{Rii* R$i\ 



(28b) 



With the fl-matrix algorithm, if both u<" +1 > and are 
needed, the full matrix recursion has to be performed even 
when u<°> = 0. If a 101 is not needed and ^ = 0, then 
quarter B-matrix recursion with Ru* is possible, but the 
R matrix has to be initialized by 



*--[! -J 



* (29b) 



With this initialization, R% p) =* Rn « 0 and R$ « -1 
for all p. 

Finally, we shall provide operation counts per grat- 
ing layer for the variants of the two recursive matrix al- 
gorithms that have been presented in this paper. The 
operation counts will be given in units of flops. 21 The 
! counts do hot include the effort in assembling the W ma- 
, trices and in solving the final linear system to yield 
• and d tW . For convenience, we shall consider only the 
' operations that are proportional to iV 3 , where N is the 
1 truncation order, the dimension of the submatrices. The 
method of counting is based on well-established rules 2 *: 
suppose that A, B. and C are N x N nonsparse matri- 
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Algorithm Number 

la 
2a 
3a 
4a 
lb 
2b 
3b 
4b 
5b 



Table 1. Operation Counts (in JV 3 - 
Different Variants of the S-Matrix 



Algorithm 



Stability 



W -o $ — > 5 
W -*> t S 
)V -> i -> r r ^ J5 
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Flops) per Grating Layer for 
and JMVtatrix Algorithms 



Operation Counts 



Unconditional 
Unconditional 
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ces; then AB + C. A'\ and A^B + C take N 2 , N 3 , and 
(4/3)N 3 flops, respectively. 

The results are summarized in Table 1 where imple- 
mentation variants of the S-matrix and i*-matrix algo- 
rithms that have been described in Sections 3, 4, and 5 are • 
represented symbolically. For example, W t ->S 
represents the variant of the S-matrix algorithm that 
uses Eq. (7) to compute the interface t matrix and then 
uses Bqs. (19a) to perform the S-matrix recursion. The 
broken arrows indicate that in some grating models the 
t matrices are obtained without using the W matrices. 
In this case, (32/3W* flops should be subtracted from 
he operation counts in Table !. The subheadings of 
the last three columns stand for full-, half-, and quarter- \ 
matrix recursions, respectively. Since the half-matrix . 
recursion of the R matrices serves no useful purpose, the 1 
corresponding operation counts are not given. Clearly, 
algorithms 2a and 3a are the most efficient, assuming 
that we start with the W matrices. 

7. REMARKS 

A. Algorithm of Chateau and Hugonin 
It is easy to see that the algorithm proposed by Chateau 
and Hugonin, 8 except for the notational differences, is 
algorithm 2a in Table 1 for the special case in whichu - ; 
0. It is one of the most efficient variants of the general 
S-matrix algorithm, but it can be slightly improved. In , 
Ref 8 each of the three factors of the layer t matrix 
[see Eqs. (6) and (7)] is passed through the recursion, 
formula separately. So the operation count, including 
the inversion ofW^, is (50/3)N 3 for the 1«*J*^*- 
In comparison, the use of product t = W 
Eq. (19a) costs 15AT 3 flops. 

B. B-Matrix Algorithm and Differential Method 

For the differential method it is natural to use the 
B-matrix algorithm because here the V -V basis is the 
natural basis. The factorization of the layer t matrices 
is unavailable in the differential method, so the applica- 
tion of the ^-matrix starts with the layer t matrices. As 
explained in Sections 4 and 5, use of the i matrix in the 
^matrix algorithm makes the stability of the algorithm 
conditional. Although the modal methods were assumed 
when we analyzed the cause of numerical instability 
of Eq, {Mb*), the conclusion applies to the differential 
method as well, because the basis functions are asymp- 



totically the same in the two cases. The first few rows of 
Tables 1, 2, and 3 in Ref. 9 clearly indicate that if the nu- 
merical layer thicknesses are not kept low, the J?-matrix 
algorithm fails when applied to the differential method. 

The R -matrix algorithm that is used in Ref , 9 is 
algorithm 4b in Table 1 of this paper, which takes 
(3lfc)N* flops per layer for the full-matrix recursion. 
lt.can.be improved slightly by using algorithm 5b, which 
takes (25/3)N a flops per layer. Instead, if the u-d ba- 
sis [functions and algorithm 4a are used, significant im- 
provement can be achieved. The operation count for the 
S-matrix algorithm is only (13/3)7^ flops per layer for 
thelhalf-matrix recursion. Furthermore, the extra work 
of solving the final system of linear equations, Eq. (17b), 
is avoided. 

C. Ji-Matrix Algorithm and Classical Modal Method 
In the classical modal method, 7 thanks to the orthogonal- 
ity [of the modal functions and the fact that the pivotal 
submatrix 4? ~ 0, not only are the t matrices obtained 
analytically without the W matrices, but the r matri- 
ces' can also be determined analytically from the t matri- 
ces' without numerical inversion of the pivotal submatrix. 
The result is the most efficient variant of the J? -matrix 
algorithm, which can be symbolized simply as r R 
Only one of the four eubmatrices of takes N 3 flops to 
construct; the rest involve only ^processes. Thus the 
overall operation counts are only (22/3)iV a and (10/3)2^ 
per layer for the full- and quarter-matrix recursions, 
respectively. 

D. S-Matrix Algorithm and Classical Modal Method 
The S-matrix algorithm can be applied to the classi- 
cal modal method, the most efficient variant being algo- 
rithm 2a of Table 1. The combination of the S-matrix 
algorithm and the classical modal method has a peculiar 
problem, which I shall describe below. 

ffn the classical modal method, the t matrices that use 
the exponential basis functions can also be obtained ana- 
lytically without the W matrices, but the s matrices in 
general cannot. Without going into any detail, suffice it 
to, say that the elements of r^ 1 at a numerical interface 
are all of the form 

(30) 



where the integration is over one grating period, ft 
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and tf/l p) are modal functions in layers p + 1 and p, {re- 
spectively! and f p is a function that depends only on the 
permittivity distribution of the two layers. In what fol- 
lows, we consider the evaluation of Bq. (30) under a spe- 
cific set of conditions: (1) the permittivity distributions 
in two adjacent layers ere symmetrical with respect to 
the origin of the x axis, (2) the grating is in the first-order 
Littrow mount/and (3) N - 2M + 1, where M is a natu- 
ral number. Under condition (1), f p is a symmetrical 
function. Under condition (2), and <ftn p) are either 
symmetrical or antisymmetries] functions, Tims if inte- 
gers I and n correspond to modal functions of different 
parities > the corresponding t matrix element is identi- 
cally zero. In the classical modal method, one normally 
indexes the eigenvalues in the order of increasing absolute 
values. Let Nc p) and ivi p) be the numbers of even and 
odd eigenvalues, respectively. Numerical experiments 
show that, under condition (2) and for a given truncal ion 
order N ~ Nl p] + Nl p \ Ni p) and Nl p} never differ, by 
more than L Thus under condition (3), either Ne P = M 
and Ni p) = M + 1, or N { * p) - M + 1 and Nl p} - M. It 
can be easily shown that if Nl p) * Ni p ^ } then all sub- 
matrices of i f P K in particular the pivotal submatrix fjj . 
are mathematically singular. Numerically, the condi- 
tion Ni p} * Nl p ~ Xi does often occur; therefore algorithm 
la of Table 1 cannot be applied to the classical modal 
method when conditions (1), (2), and (3) are met simul- 
taneously. It can be verified numerically that the (ma- 
trix sum that is to be inverted in Eqs. (19a) sometimes 



tve 



becomes numerically ill-conditioned under the abo 
three conditions; therefore algorithm 2a fails too, even 
though it does not involve the inversion of the pivotal 



submatrix f^ J - 



Lnmatnx - i 
Fortunately, the singular matrix problem described 
above can be easily avoided by the use of an even trun- 
cation order. If N = 2Af, then the characteristics of 
the eigenvalue distribution automatically guarantee jthat 
&i p) Ni p) = M. | . 

Another interesting difference between the 2Z-matnx 
algorithm and the S-matrix algorithm, as they are applied 
to the classical modal method, is that the law of enjergy 
conservation (in the case of dielectric gratings) is satisfied 



automatically by the former, but it is satisfied only 
increasing truncation orders by the latter 



with 



£. Other Possibilities j m 

The essence of the B-matrix and S-matrix algonthms is 
to avoid tha presence of the growing exponential functions 
in the matrix manipulations. In this spirit, several other 
stable algorithms have recently been presented in the 
literature. For example. Montiel and Neviere 9 presented 
an algorithm that they called the iZ'-matrix algorithm. 
In view of the current paper, it can be considered an 
S-matrix algorithm that uses the V-V basis functions. 
In Ref. 10 the i?-matrix algorithm was used with the 
exponential basis functions (maybe it can be called the 
S'-roatrix algorithm?)- The scattering-matrix approach 
of Cotter e* aL n is essentially algorithm 2a of TfeUe 1, 
except that their i matrix is the inverse of the l ma trix 
in this paper. Clearly there are many other possibilities, 
but it is pointless to enumerate all of them. j 
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8. SUMMARY 

The mathematical formulations of the S-rnatrix and 
'^-matrix algorithms have been systematically presented- 
The presentation is given in a unified fashion, indepen- 
dent of underlying grating models, grating geometries, 
and grating mountings. The physical interpretations of 
the algorithms are illustrated. In addition, many vari* 
ants of the algorithms are presented and their numerical 
stabilities and efficiencies analyzed. 
\\ The S-matrix and E-matrix algorithms are inherently 
stable because thay avoid the appearance of the expo- 
' nentially growing submatrices in the recursion formulas. 
However, to further ensure that the algorithms be uncon- 
ditionally stable, effort should be made to avoid the expo- 
nentially growing submatrices in the intermediate steps, 
i.e., in the calculation of the layer s or r submatrices. 
{Whenever the factorization, as given in Eq. (6), of the 
' layer f matrix is possible, the interface t matrix should 
' be used directly in the constructions of layer s (or r) ina- 
' trix or in the nonnormalized S-matrix (or R -matrix) re- 
•cursion. "When factorization is impossible, the S-matrix 
and J2-matrix algorithms are stable under the condition 
that the layer thicknesses and the truncation order be 
•kept low, as quantified at the end of Section 4. 
1 The comparative study of the two matrix aJgorithms 
presented here seems to favor the S-matrix algorithm. 
.The physical interpretation of the S matrix in terms of re- 
flections and transmissions is more intuitive than that o* 
the R matrix in terms of die impedance and admittance* 
.The exponential basis functions adopted by the S-matrix 
algorithm are numerically much easier to handle than the 
' trigonometrical basis functions adopted by the R matrix 
' algorithm. Based on the operation counts, the S-matrix 
algorithm is more efficient than the R -matrix algorithm. 

Many implementation variants of the algorithms are 
presented in this paper. The variants that use all in- 
termediate matrices, algorithms la and lb in Table 1, 
, are the least efficient ones. They have only pedagogical 
value. The variants that bypass some of the interme- 
' diate matrices, for example, algorithms 2a, 8a, and 2b, 
, are the most efficient ones. However, as exemplified in 
Section 7, which algorithm and variant are more efficient 
often depends on the grating model being used. It is the 
hope of the author that the information provided here will 
enable the reader to apply the most efficient algorithm to 
the grating model at his or her disposal. 

APPENDIX A 

To derive Eq. (13b), let us imagine that layer p is a sum 
of two layers, the first layer has zero thickness, with 
the layer t and r matrices given by Eqfi. (7) and (14b), 
respectively. The second layer has thickness A pi but it 
does not cross a material boundary- Its equivalent layer 
t matrix is just <f> ip K given by Eq. (4b). Denoting the 
equivalent layer r matrix corresponding to $ KP} by r • 
from Eq. (14b) we have 

( _ f-uSf 1 cotU^M nl p) csc(A^,)~l . juj. 
' " csc(A!n%) n!^ot(Ak%)J 

Equations (15b) can be viewed as a set of rules that 
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(A2) 



and matrix in relation 



(A3) 



to obtain matrix fl" 0 in relation 



Here we have . i 

i 

[ [/<">(* -o) j _r J ! 

and what we want is the matrix in 

Through comparison of Eq B . (A2)-CA4) with Eqs, (A6)- 
(A7), it is evident that the two sets of equations have 
identical algebraic structures. Therefore Eqs..(18b) can 
be obtained from Eqs. (15b) provided that H», r», and 
?W are identified with f<», Jl««- I », and fr*, respectively. 



APPENDIX B 

When t\£ = 0, which happens in the classical modal; 
method, Eq. (14b) cannot be used to compute but ny 
this case for sure r&'^O. Suppose that r^' is noftsinguj 
lar and aS»" h p + W, then Eq. (14b') can be used to derive 
the layer r matrix. After some simple algebra, we havej 



c A,)~|i 

X) Jl 



["[# " W cottA^)]^- 1 «W» 
L CBli) 

■ i: 

This expression off"", like Eqs. (13b) contains w .expo- 
nentially growing functions, and therefore is suitable fur 
the unconditionally stable JJ-matrix recursion. 

APPENDIX C i 
Similarly to the treatment in Appendix A, we consider 
that each of the two factors in Eq. (6) corresponds to ja 
layer, one with zero thickness and the other with the fun 
tMekness/u. Substituting forPrtinEois.(20b),ana 
making some simple algebraic rearrangement, we obtajn 



an intermediate J? matrix, which is given by 
I'kf - - Hif cot(A«{"A p ) + .„<>> cscU\>%) 



where 



w (p)=[^cot(A^A p )-2?{rV, 



(CD 



(CI') 



Clearly, all eubmatrices of R ip) are of order Oil) and 
numerically stable. To complete the nonnormalized 
jR -matrix recursion using the interface t matrix, we need 
onlyito use Eqs. (20b) again, this time replacing R ip ~ l) 
by fi* p) and by t lp) . The result is 
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